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ABSTRACT 


The goal of this research was to determine kinematic parameters of the lower limbs of a 
subject pedaling a bicycle. An existing measurement system was used as the basis to develop 
the model to determine position and acceleration of the limbs. The system consists of an 
ergometer instrumented to provide position of the pedal (foot), accelerometers to be attached to 
the lower limbs to measure accelerations, a recorder used for filtering, and a computer 
instrumented with an A/D board and a decoder board. The system is designed to read and 
record data from accelerometers and encoders. Software has been developed for data collection, 
analysis and presentation. Based on the measurement system, a two dimensional analytical 
model has been developed to determine configuration (position, orientation) and kinematics 
(velocities, accelerations). 

The model has been implemented in software and verified by simulation. An error analysis 
to determine the system’s accuracy shows that the expected error is well within the specifications 
of practical applications. When the physical hardware is completed, NASA researchers hope 
to use the system developed to determine forces exerted by muscles and forces at articulations. 
This data will be useful in the development of countermeasures to minimize bone loss 
experienced by astronauts in micro gravity conditions. 
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INTRODUCTION 


Decrease in muscle loading and external loading of bone during weightlessness in space 
can result in cancellous bone loss of 1% per month in the lower extremities and 2% per 
month in the calcaneous (Figueroa, 1993). Therefore bone loss is a serious problem encountered 
by astronauts that must remain in micro gravity conditions for the duration of a mission. It is 
hypothesized that loading bone appropriately during exercise may prevent bone loss (Figueroa, 
1993). To minimize bone loss, NASA scientists are considering development of exercise 
countermeasures. This development involves definition of exercises and doses that will stress 
the bone so as to minimize bone loss on subjects participating in a bedrest study. The bedrest 
model is used to simulate micro gravity conditions. Loads on exercises that are considered 
effective must be quantized to be compared with loads measured during exercise in space. Thus, 
information about the kinematics, dynamics and forces exerted by a particular muscle or 
muscle group during exercise is necessary to quantify bone loading for the development of 
exercise countermeasures. 

This thesis describes a system to determine the kinematic parameters of the lower limbs 
during exercise. The system elements were selected based upon specification from NASA. It 
consists of instrumentation, sensors and methodologies necessary to support the determination 
of kinematic and dynamic information. It encompasses a combination of small, light, and robust 
sensors, suitable for use in the confined environment of a space vehicle. 

The research included refurbishing and improvement of existing hardware, development 
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of software, and development of a two dimensional model to determine kinematic parameters 
of the lower limbs of a subject riding an ergometer. 

The equipment used includes an ergometer, a 28-channel recorder (XR-9000 Cassette 
Data Recorder, Teac Corp., Tokyo, Janpan), accelerometers (EGAXT-10, Entran Devices Inc., 
Fairfield, N.J.), three rotary encoders, a 486 class personal computer fitted with two cards: (1) 
a three-channel decoder card (Model 5312 board, Technology Inc., Minneapolis, MN) to read 
the encoder information and (2) a 64 channel A/D board (AT-MIO-64F-5, National Inc., Austin, 
Texas) to digitize the data from the accelerometers. 

Accelerometers are attached to the calf and thigh sections during exercise. Signals from 
the accelerometers are amplified before recording. Filtering is done using the recorder. The 
angular positions of the crank and the pedals may be measured by encoders (M2005 122 1031, 
Dynapar Corp., Gurnee, IL) installed at the joints of the crank and pedals of the ergometer. 
The recorder performs anti-aliasing filtering. Figure 1.1 shows the hardware connections. 

Data from the accelerometers and encoders are synchronously recorded. Data from 



Figure 1.1 






3 


encoders is converted to angular positions, angular velocities, and angular accelerations. Data 
from the accelerometers is used in the kinematic analysis. The A/D board used in the research 
is the AT-MIO-64F-5 board produced by National Instruments Inc. It is a multifunction analog, 
digital, and timing I/O board for the personal computer. The board has a 5 /xsec, 12-bit 
sampling ADC that can monitor a single input channel, or scan through the 64 single-ended or 
32 differential channels. NI-DAQ application programming software was used to drive the AT- 
MIO-64F-5 board. 

A novel two dimensional analytical model was developed, which uses data frorn the 
above system as input to determine kinematic parameters. The model has been implemented in 
a computer program, KINET.M. The kinematic parameters are needed to determine bone and 
muscle loading, but this thesis does not include determination of loads. 

Section Two describes the theoretical background of this research. Section Three gives 
a complete description of the theoretical model for kinematics. The error of the angular positions 
of the lower limbs resulting from the precision errors of the accelerometer values is also 
analyzed. Section Four describes the software of the proposed system. The last two sections 


include conclusions and recommendations. 



BACKGROUND 


According to Wolff ’s law, living bone changes according to stress and strain conditions 
on the bone (Fung, 1919). The mechanism of bone remodeling is linked to strain or change in 
bone dimensions caused by stress applied to it. To increase bone strength, stress must be 
applied to produce sufficient strain. To develop exercise countermeasures that produce bone 
stress adequate to minimize bone loss in micro gravity conditions, it is necessary to quantify the 
loads during exercise. 

Generally, forces and torques exerted at joints are measured by use of the equations of 
motion (Newton-Euler Law of Motion) of the body parts of interest (Figueroa, 1993). This 
results in a set of equations that relate the torques and forces applied by muscles, and by 
neighboring bone sections at the joints (Redfield and Hull, 1986; Anderson et. al., 1993; Yang 
et. al., 1993; Abdel-Rahman and Samir, 1993; Verstraete, 1991; Ericson et. al., 1985; Harrison 
et. al., 1986). To use the equations of motion, the measurement of kinematic parameters such 
as accelerations, positions, and orientations of the body parts is required. This is the topic 
addressed by this thesis. 

Popular methods to measure the kinematic parameters during exercise use camera/light 
systems (Ericson et. al., 1985; Harrison et. al., 1986). These methods are accurate, but the 
instrumentation is bulky and requires large spaces, and the procedures to install and/or calibrate 
are time consuming. The systems also require unobstructed line of sight between cameras and 
lights, and are generally expensive. Another method, to describe the kinematics of the leg in 
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the sagittal plane only, assumes that the hip remains static, and applies kinematics of a five-bar 
linkage (Redfield and Hull, 1986A, Redfield and Hull, 1986B). The system developed uses 
classical engineering methods, more in line with the requirements put forth by NASA (compact, 
transferable to a space environment, easy to install and operate), and uses accelerometers and 
position sensing instruments. The method described in this thesis is unique in that although it 
uses accelerometers, position is not obtained by integration. This avoids the usual accumulation 
of errors characteristic of methods that integrate measured accelerations to determine position. 
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THE PROPOSED SYSTEM 


3.1 INTRODUCTION 

The proposed system uses a combination of accelerometers and rotary digital encoders. 

After passing through anti-aliasing filters, measurements made by these sensors are read into a 

computer with data acquisition cards: (1) AT-MIO-64F-5, National Inc., Austin, Texas and (2) 

Model 5312 board, Technology Inc., Minneapolis, MN; and then those measurements are used 

■* ^ 

to determine the configuration and kinematics. 

Section 3.2, Methods, and Section 3.3, Kinematic Equations, describe the measurement 
system and the methodology. Section 3.4 describes an analysis of error propagation from the 
accelerometer values to the angular position of the body section. 

3.2 METHODS 

To apply the Newton-Euler Equations of Motion, the acceleration of the center of mass 
must be first determined. To obtain acceleration of the center of mass, miniature accelerometers 
are used to measure the acceleration vector of two points (for motion in the sagittal plane). The 
acceleration vector will be expressed with respect to a coordinate system attached to the 
member. The orientation of this coordinate system will be determined using the acceleration of 
a known point in both coordinate systems (a base coordinate system or the inertial frame, and 
the one attached to the member). The relative rotation of the two coordinate systems will be 
determined by the equations that relate the orthogonal components of the same acceleration 
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vector expressed in the two coordinated systems rotated with respect to each other. 

Using accelerometers on Earth requires that the component of the acceleration of gravity 
along the direction of the accelerometer axis should be included as part of the measured 
acceleration. This component will be considered in the development of the model described in 
the next section. 

3.3 KINEMATIC EQUATIONS 

The formulation assumes that the acceleration of one point of the member of interest 

r * **• 

is known, and that accelerometers placed on the member measure acceleration components along 
the directions of the frame attached to the member. For example. Figure 3.1 shows the thigh, 
calf, and foot projections onto the sagittal plane. The position and orientation of the foot are 
measured using optical rotary encoders attached to the crank and pedal rotation axes. Angular 
velocities and angular accelerations of the crank and pedals can thus be known. 



Figure 3. 1 
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3.3.1 Determination of angular velocity and acceleration 

The acceleration of the foot joint with respect to the inertial frame, denoted °a 2 , can be 
determined using Equation 3.1. In this equation l { is the length of the crank, l 2 is the length 
of the pedal/foot, 0, is the angular position of the crank w.r.t the inertial frame, 0 2 is the relative 
angular position of the pedal with respect to the crank joint, w, and a, are the angular velocity 
and acceleration of the crank, w 2 and a 2 are the relative angular velocity and acceleration of 
the foot with respect to the crank, Cj is cos^), S, is sin(0[), C 12 is cos(0j+0 2 ) and so on. 


V 


/jCi ^2^12 


«1 

2 2 
<Xi+a 2 


k^i2 
i ^12 


(<^>i) 2 

(0) j +0) 2 ) 2 


(3.1) 


Two accelerometers are attached to the calf; and another two to the thigh. The 
accelerations measured by these accelerometers are used to determine the angular velocities, 
angular accelerations, and linear accelerations of the centers of gravity of these sections. Since 
the accelerometers are used on Earth, a lg value in the vertical direction is added to the actual 
accelerations. As shown below, the gravity component increases the complexity of the method 


for determining the angular positions of the calf and thigh sections. 

The linear acceleration of the knee joint with respect to Frame 3 is denoted 3 a 3 . The 
linear acceleration of accelerometer 1 attached to the calf with respect to Frame 3, is denoted 
3 a cl , and can be expressed as follows: 


\r 3 a 3 +a 3 x 


(3.2) 


where a 3 denotes angular acceleration of the calf, to 3 angular velocity of the calf, and r cl 



Figure 3.2 


position vector of accelerometer 1 on the calf with respect to Frame 3. 

Because all motions are assumed to be on the sagittal plane. Equation 3.2 can be written 
as follows in two dimensions. 


3 r „ _ 3- 2 _ 3 _ 3 - (3.3) 

r cly tt 3 r clx“> 3 “ a clx “ax V 


3- ~ _ 3_ ., 2 _ 3 fl _ 3^- (3.4) 

r clx a 3 r cly° ) 3 a cly °3y 


Equation 3.3 and Equation 3.4 can be expressed in matrix form as Equation 3.5. 
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\,=V i 


r ci> r cijc 


r , - - r 

clx 


cly 


a, 


0)< 


(3.5) 


Similarly, the acceleration of accelerometer 2 attached to the calf, denoted 3 a c2 , can be 
expressed as follows: 


3 r^,a, - 3 r.^co, 2 = 3 a c2x - 3 fl 3x 


, c2y u '3 c2x 3 


(3.6) 


3 T.,jx, - 3 r . o) 2 = 3 a r2v - 3 a 3 


(3.7) 


r c2* u 3 '^3 _ v *c2y ^ 

Equation 3.6 and Equation 3.7 can also be expressed in matrix form as Equation 3.8. 


3fl c2=Vi 


- 3 r 


c2y r clx 


clx 


'c2yj L' 


a. 


0), 


(3.8) 


Before evaluating cj 3 and o; 3 , the relation between the accelerometer readings and the true 
accelerations must be established. As shown in Figure 3.2, suppose that the rotation angle 
between Frame 3 and the inertial frame is 0 O3 - Gravity exerts an acceleration of -gsin( 0 O3 ) in 
the direction of the x-axis of Frame 3 and an acceleration of -gcos( 6 0 3 ) in the direction of the 
y-axis of Frame 3. Therefore, the relation between the reading of accelerometer 1 attached to 
the calf, denoted a cl , and the actual acceleration, denoted 3 a cl , is as follows: 

\ Jx = a ci* + Sswflra) (3 ' 9) 

\iy = a dy + SC 0 S(Q 03) (3 - l0) 


or in matrix form 
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|#sm(0 O3 )l (3.11) 

3 o ,=o ,+ 

C C [gCOS(0 O3 ) 

Similarly, the relation between the actual acceleration of the point of attachment of 
accelerometer 2, denoted 3 a c2 , and the reading of accelerometer 2 attached to the calf a c2 is as 
follows: 

= a cix + 8sin(Q 03 ) ( 3 - 12 ) 

\ly = a c2y + 8COS(Q 03 ) ,J3.13) 

or in matrix form. 


\ 2 =a c2 + 


gsin(d Q3 ) 

gcos(d 03 ) 


Until the rotation angles between the inertial frame and Frame 3 or Frame 4 can be 
determined, the true accelerations of the accelerometers attached to the calf or thigh sections can 
not be obtained. To simplify the equations, a subscript g will be added to denote an acceleration 
to which a corresponding gravity component has been added. For example, 3 a 2g = 3 a 2 - J g. 

Substituting Equation 3.9 and Equation 3.10 into Equation 3.3 and Equation 3.4; 
Equation 3.12 and Equation 3.13 into Equation 3.6 and Equation 3.7, respectively. The newly 
defined equations are: 


3 r rN a 3 - 3 r . co 3 2 = a, + gsin(Q 03 ) - V (3-15) 
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\l x « 3 - \ly<* 3 2 = + ^«(e 03 ) - 3fl 3j (3 ’ 16) 

" " 3r cW = a c2* + Ssin(Q 03 ) - ^ (3 - 17) 

3r c2x a 3 - \2y“> s' = «c 2 y + ^( 603 ) " S, (3J8) 

Subtracting Equation 3.15 from Equation 3.17 and Equation 3.16 from Equation 3.18, 
respectively, and expressing the results in matrix form, 

-(\ 2 y-\iy) ~( 3r c2x~\J 
(\2x-\lx) -(\ 2 y-\ly\ 


a. 


co, 


a c2 a cl 


(3.19) 


Because 3 r cl and 3 r c2 are known anda cl , a c2 are measured by the accelerometers, 
angular velocity and angular acceleration can be evaluated as follows: 


-(\2y-\ly) ~ (3r c2x^ r clx) 
(\2x-\lx) <\2y~\ly). 


( a c2~ a cl) 


a. 


to, 


(3.20) 


3.3.2 Determination of Angular Position 

Once a 3 and co 3 are evaluated, Equation 3.15 and Equation 3.16 can be rewritten as 


follows: 
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a clx <- \ ly a 3 - \i^ 3 2 )= K- ^(0 O3 ) = 3 «3^ (3 21) 

a cly< \li a 3 - 3r cly“>3 2 )= ^( 6 03 ) = ^ ( 3 - 22 ) 

All the terms on the left-hand side are known; therefore, the terms on the right-hand 
side, denoted 3 a 3g -meaning 3 a 3 with the added gravity component- can be determined. 

From rigid body kinematics, the relation between the accelerations of the knee joint J a 3 
and of the foot joint 3 a 2 is given by: 


3 r j. 2 + 3 _ _ 3_ (3.23) 

r 32y ft 3 r 32i W 3 a 3»- °Zt 


3*. ~ _ 3_ (i) ^ (3.24) 

'32i a 3 r 32y W 3 a 3;y a 2;y 

where 3 r 32 stands for the position vector of the foot joint with respect to Frame 3. 

Adding -gsin( 0 O3 ) to both sides of Equation 3.23 and -gcos( 0 O3 ) to both sides of 
Equation 3.24, respectively, one obtains, 


- 3 r 32y a 3 - V 32 x co 3 2 + - gsin(Q 03 ) = z a^- gsin(6 0 3 ) (3-25) 

3r 32* a 3 ‘ 3r 32y t0 3 2 + \ ~ SCOS(6 Q3 ) = \ " gCOS( 6 Q3 ) (3-26) 


If 3 a 2 plus its respective gravity component is denoted as 3 a 2g , the above two equations 


become: 
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3r 32v tt 3 r ^7x S) ^ + fl 3gx a 2 j 


■gx 


2.3, 


V„ v a, - + "a 3 ^- '’a 2 «y 


(3.27) 

(3.28) 


32x“3 ' 32y 3 

Since the terms on the left-hand side are all known, the ones on the right-hand side, "a 2g , 
can be determined. 

From robot kinematics theory (Phillip John McKerrow, 1991), the relation between the 
acceleration of the foot joint with respect to the inertial frame and the acceleration with respect 
to Frame 3 is as follows: 


"a* 





sin^) cos(0 O3 ) 


a 


2x 


2yJ 


(3.29) 


If 3 a 2 is expressed in terms of 3 a 2g , then Equation 3.29 becomes Equation 3.30. 




. V 



sin(0Q 3 ) 
sin(0 O3 ) cos(0 O3 ) 


3 a 2gx +gsin(Q 03 ) 

3 a 2gy +gcos(Q 03 ) 


(3.30) 


In Equation 3.30, the only unknown variable is 0 O3 , and it can be expanded as 


follows: 


°a 2x = 3 a 2g3 .cos(0 o3 ) + gcos(0 O3 )sin(0 O3 ) - 3 o 2gy sin(0 O3 ) gcas(0 O3 )sm(0 O3 ) (3. 31) 


°a 2y = 3 a 2 ^sin(0 O3 ) + gsin 2 (d 03 ) + 3 a 2g> cos(0 O3 ) + gcos 2 (6 03 ) 


(3.32) 


Then the two equations become: 



(3.33) 


V = 3 o 2 ^cos( 6 03 ) - 3 a 2 ^sin(0 o3 ) 

\ y = 3 a 2 ^sin(0 O3 ) + 3 a 2 ^cos(0 O3 ) + g (3 ' 34) 


Finally d 03 can be evaluated as: 


i u 0 y 

0 O3 = cos- [ 


3fl 2 gx + a 2y a 2gy 8 tt 2gy j 

(V 2 + ( 3 V 2 


(3.35) 


3.3.3 Determination of Linear Acceleration 

The true acceleration of the knee joint with respect to the inertial frame can thus be 
evaluated in two steps. First evaluate the true acceleration with Fespect to Frame 3: 


X = 3 <V + ^* (0 O3) 

\ = 3fl 3 o + SCOS(Q 03 ) 

Then the true acceleration of the knee joint with respect 
by ao 3 , is obtained as follows: 


(3.36) 

(3.37) 

to the inertial frame, denoted 


V 


cos(0 O3 ) -sin(0 O3 ) 



° a 3y 


sin(0 O3 ) cos(0 O3 ) 




(3.38) 


After finding a^, the same procedure is applied to determine the kinematic parameters 


of the thigh section using the acceleration of the knee joint. 
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3.4 ERROR ANALYSIS 
3.4.1 Analysis of error propagation 

Accelerometers produce errors in precision and accuracy. If the accuracy error is not 
negligible, then the accelerometer must be recalibrated. It is assumed that errors are due to 
precision alone, so that the resultant error of angular position of the moving frame attached to 
the body member is caused by the propagation of the precision error of the accelerometers. An 
analysis of accuracy error propagation can be inferred from the precision error analysis and will 
not be considered in this thesis. 

The uncertainty of the measured accelerations is denoted by 5. The error analysis follows 
the procedures used to determine 0 O3 and 

The only error carrying inputs to determine d 03 in Equation 3.39 are 3 a 2g and °a 2 . The 
latter is obtained through calculations using encoder values. The encoders (M2005 122 1031, 
Dynapar Corp., Gurnee, IL) are of good resolution (2048/rev). It will be assumed that the only 
sources of error of 6 03 are the values from the accelerometers, a clx , a cly , a c2x and a c2y . 

s 3 a 

° (3.39) 

Define k as follows: 

°^2x 3 < 3 2 gx + °% 3 fl 2 gy ~ 8 a 2gy (3 40) 

( 3 v 2 + 

where 3 a 2g is a function of a 3 , co 3 2 and 3 a 3 as follows, 


e 03 = cos ' 1 


V + \ a 2gy 


(\f + 



( 3 . 41 ) 


3fl 2 gx 3r 32y a 3 r 32x^3 + U 3gx 


3/i —3 r _ 3- y \ 2 3« (3.42) 

r 32x a 3 r 32y (j *3 a 3gy 


where 3 a 3 is also a function of a cl , a 3 and o> 3 2 , 


3/1 — n —(— 3 r cc — 3~ q 2 \ (3.43) 

a 3#x“ V r c;^ a 3 r c^ w 3 J 


V' a *,< ~ ) (3 ' 44) 


and a 3 and cj 3 2 are both functions of a cU , a cly , a c2x and a c2y . 


a. 


r c21y( a c2x a cb) + 3/ c2Jj( a c2y g cip 


3,2 


3,2 


^c 2 ijc + ^c 2 i;y 


(3.45) 


a). 


r c2lx^c2x a clx ) r c21y^c2y &cly) 


3 2 
.j r 


3,2 


c2 lx + ^"c2 ly 


(3.46) 


Thus 3 a 2 , x and 3 a 2gy can be expressed in terms of a clx , a cly , a c2x and a c2y as follows. 


3 a 2gx ={ [ ( 3 r C 2ix 2 + 3 r C 2i y 2 )+ 3 r c 21 y ( 3 r cly - 3 r 32 y )+ 3 r c 21 x ( 3 r clx - 3 r 32x ) ] a c)x + 

[ 3 r c 2 i y ( 3 r clx - 3 r 32 x )- 3 r c 21 x ( 3 r cly - 3 r 32y ) ] a cly + 

[- 3 r c 21 y ( 3r cl y - 3 r 3 2 y )- 3r c21x( 3r clx- 3r 32x) ] 3 c2x + 

[ 3 r c2 ,x( 3 rcl y - 3 r 3 2 y )- 3r c21 y ( 3r clx- 3r 32x) 1 a c2y 

} / [ 3 r c 2 ix 2 + 3 r c2i y 2 ] 


(3.47) 
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3 a 2gy ={ [- 3 r c 21 y ( 3 r clx - 3 r 32 x )+ 3 r c 21 x ( 3 r cly - 3 r3 2y ) ] a clx + 

[ ( 3 r c 2 ix 2 + 3 r c 21 y 2 )+ 3 r c 21 y ( 3 r oly - 3 r 32 y )+ 3 r c 21 x ( 3 r clx - 3 r 32x ) ] a cly + 

[- 3 rc 2 ix( 3 r c i y - 3 r 32y )+ 3 r c21y ( 3 r clx - 3 r 32x ) ] a c2x + 

[* 3 r c 2 i y ( r oly - r 32y )- r c2 [ x ( r clx - r 32x ) ] a c2y 
} / [ 3 r tI1I 2 + 3 r c2ly 2 ] (3.48) 

where 3 r c21 = 3 r c2 - 3 r cl is the vector from accelerometer 1 to accelerometer 2 expressed in Frame 
3. 


0 O 3 is defined by: 


0 O3 = cos~ l (k) = cos _1 [ k( 2 a 2gx , (3-49) 

which may be expressed as follows: 


0 O3 =cos _1 ( k[ 


^Igx^clx^cly^clx^cZy ^ * ^Jgy^clx^cly^clx'^cly) 1 ) 


(3.50) 


The square of the uncertainty of 0 O3 is defined by: 



fa chc da cly da c2x da c2y 


(3.51) 


Using the chain rule, the partial derivative terms of the preceding equation are 
determined. 


dk d Q-2gx + ^03 dk d 
da clx dk d 3 a 2gx da clx dk d 3 a 2gy da cIx 


(3.52) 
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^®03 __ 

3003 0* 

3 3 ^2gar + 3003 3& 


^ a cly 

3* a 3 a 2 ^ 

0a c7 , 3fc a 3 a 2 ^ 

3a c/^ 

00 O3 _ 

^®03 dk 

3 \ x H 30 O3 a* 

e3a *a 

^ a c2x 

dk 3 3 0 2 ^ 

da cix dk 3 3 a 2 ^ 


Q_> 

CD 

o 

u> 

II 

30q3 dk 

3 3 <V , 30 O3 0* 

a \» 


da c2y dk a \ gx da c2y dk a z a 7gf da c2y 


(3.53) 


(3.54) 


(3.55) 


where 


00 eg 

~dk 


1_ 

\[i - k 2 


(3.56) 


dk 


0fl Zc( a 2gx + °2gy) ^ a 2gx + *hy ^2gy & 


0 3 a 




( 3 4* +3 <4) 2 


0 3 


0£ _ ( ° a 2y ^)( 3a 2gJ +3fl 2gy) £ 3a 2gy(° a 2i ^hgx^fhy 3g 2gy # fl 2gy) 

( 3 4* +3fl i) 2 


(3.57) 


(3.58) 
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2 S*_ 


0 a 


clx 


( + 3 4/>) 

+ 3r c21y^ r cly r 32p + r c2// r cix 


3r c2/x + 3r c2ij 

a3 ^_ 3 W 3 ' 

r«i,- ‘fJ-UV ’W 

d a cly 

3-2 . 3_2 

r c21x r c2Jy 

d \gx_ W 3 

r „,~ ^KvAA/vT 3f 32,) 

a c2x 

" 3 2 3 2 

^c27j: + '"c2i_v 

3 W 3 

'•civ" W -\ 2 ly( 3r clx- \tx) 

3 fl c2> 

3 2 3 2 

^c21x + r c22v 

3 3 <V_ - W 

VclT 3 '-32x) +3/ -c2i,( '’civ" r 32y) 

3 fl ci;c 

3 2 3 2 

r c2ij + r c2iji 


( 3 . 60 ) 


( 3 . 62 ) 


( 3 . 63 ) 


^ 3fl 2ey _ ( 3/ "c27^ + ^ r c21y ) + ^ r c 2l y^ T cly *32 y) + K c2lx^ T clx r 32j) 


0 a 


2gy _ 

c ^y 


3 2 3 2 

' c2ij: + %2i. 


( 3 . 59 ) 


( 3 . 64 ) 
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3 3 « 2g y - \2lAcly- 3 W +3r c2/A/,- r 32x) 


d a 


c2x 


3 2 3 2 

^c2ix + ^c2Iy 


(3.65) 


3 3fl 2gy _ f c 2ly( r cly r 32p r c2ix( r c2x r 32x) 


a 


'c2y 


3 2 3 2 

^c2ix + ^c21y 


(3.66) 


In Equation 3.56, when k approaches one, the partial derivative of 0 O3 approaches 
negative infinity. This presents a practical problem to the error analysis, so the error of k is 
analyzed instead of the error of 0 O3 . 


The square of uncertainty of the k is as follow: 


2 =( « 2+( « )2(S p 2 +( * ) 2(5 )2+( * ^(5 ,2 (3.67) 

da clx da cly da c2x da c2y 

Using the chain rule, the partial derivatives terms in Equation 3.67 are defined as 


follows. 


dk dk d \ gx dk d\y (36g) 
da a x a da cu a da dx 

dk_ = dk a t dk a 3 o 2gy 
da dy a *a 2gx da cly d s a 2gy da cIy 


(3.69) 
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dk dk & 3fl 2gy + dfc d °2gy 

a °c2x 3 3 a 2gx da c2x d 3 a 2gy da c2x 

dk _ dk d 3 ^ g-t: ^ dk d 

5o c 2, a 3 o 2 ^ a« c2y a 3 o 2gy ^ 

The partial derivative terms of Equations 3.68-3.71 have been determined as Equations 
3.56-3.66. 

The square of uncertainty of k is thus determined using Equations 3.67-3.71. Finally, 
the square of the uncertainties of 0 O3 can also be determined using the equation below. 


(3.70) 


(3.71) 



, ae 03 

= s- 1 (-£) 

dk 


(3.72) 


In Equation 3.72, the value of the partial derivative term affects the outcome of the 
square of uncertainty in 0 O3 . The partial derivative term has been defined as Equation 3.56. 
Within 90 degrees, it is seen that since cos( 0 O3 ) is decreasing, the square of the partial 
derivative term is also decreasing. This means that less 0 O3 results in larger value of the square 
of the partial derivative term, and thus in larger error. Therefore, cos( 0 O3 ) is not suitable for 
error analysis 90 degrees, and sin( 0 O3 ) should be used. Conversely, cos(0 o3 ) is more suitable 
than sin( 0 O3 ) within 90 up to 180 degrees. 

From Equations 3.53 and 3.54, 0 O3 can also be evaluated as: 


e 03 = sin - 1 [- 


( S - 8? Cl, 




S* Sty 


+ ( 3 V 2 


] 


(3.73) 
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Also, k should be defined as sin( 0 O3 ) and thus Equations 3.57 and 3.58 become: 

dk ( \ - *)( 3 <&+ 3 4j- 2 Vl 3 V (3,74) 

sK ( 3 4. +3< 4) 2 

dk _ “ 0a 2x( 3a 2y +3a 2gy^~^ 3fl 2 g vt 0fl 2y 3fl 2g):~ 0a 2, 3(I 2gt # (3.75) 

( 3 4 »* 3 < 4) 2 

Yet Equations 3.59-3.66 remain the same. 


3.4.2 Numerical demonstration 

Based on the analysis in the previous section, a numerical estimation of the error ot k 
was performed. The numerical data sets used as inputs are designed based on assumptions that 
angular velocities of all the sections are constant and angular accelerations are all zero. Other 
assumptions to simplify the creation of the numerical data are that the hip remains static, the 
rider-bicycle system is modeled as a five-bar linkage, and the pedal positions are always 
horizontal. The five-bar linkage system and position vectors on the calf section are shown in 
Figure 3.3. If the angular position of the crank is given, the angular positions of the other 
moving sections can be determined. At the same time, the angular velocities of the calf and 
thigh are chosen to be respectively one fourth and one sixth of the angular velocity of the crank 
chosen to be 6 rad/sec. Thus the angular velocities of the calf and thigh are respectively 1 .5 and 
1 rad/sec. Other numerical data are defined as follows: accelerometer precision errors are 
assumed to be 1%, 3%, 5% and 10%; and 3 r cljt = -0.10 m, 3 r cly = -0.05 m, 3 r c21y = 0 m. 
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3 r3?<= _ 0 .40 m, 3 r 32y = 0 m and 3 r c21y = 0 m. To see the effects of both the distance between 

the two accelerometers and the precision errors of accelerometer, 3 r c21x is chosen to be -0.20 
m and -0.3 m. Using the preceding position vectors, Equations 3.59-3.66 can be determined. 



In this numerical demonstration, two MATLAB programs TEST.M and KINET.M are 
used and their descriptions are respectively in Sections 4.3 and 4.4. The TEST.M program 
produces accelerometer values for simulating the kinematics model. The KINET.M determines 
the kinematic parameters which have been described in Section 3.3, Kinematic equations. From 
the geometry of the five-bar linkage model, is given and 6 2 is determined from the condition 
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that the foot remains horizontal, thus 0 O3 and 0 M can be determined. Because 0 O3 in this five-bar 
linkage model is less than 90 degrees, k is defined as sin( 0 O3 ). The accelerometer values are 
produced by TEST.M using the above 0„0 2 , 0 O3 , d M and the designed angular velocities of the 
moving sections. Then KINET.M determines 3 a 2g and °a 2 . Once 3 a 2g and °a 2 are known, 
Equation 3.57 and 3.58 can be evaluated. Finally, the uncertainty of k is obtained using the 
variations of the accelerometer values ( Equation 3.76). The resultant variations of 0 O3 can be 
evaluated using Equation 3.72 and the results are shown in Figures 3. 4-3.5. 

These two figures show error calculations at different distances between accelerometers. 
Each of the four curves in each figure represents different accelerometer precision errors of 1 %, 
3%, 5% and 10%. The results shown in Figures 3.4 and 3.5 are at the same crank’s angular 
velocity-6 rad/sec, but at different distances between accelerometers. 

From Equations 3.59-3.66, it can be seen that increasing the separation between 
accelerometers decreases the error. This is proven by comparing Figure 3.4 and 3.5. The 
results show that within acceptable accelerometer precision errors (<5%), the expected error 
is probably within the specifications of practical applications. 



Angular error (degrees) 




Angular error (degrees) 


Distance between the two accelerometers = 0.3m 



35 40 45 50 55 60 65 

Absolute rotation angle of the calf (degrees) 


Accelerometer precision error 10% ***** 

5% xxxxx 
3 % + + + + + 
1% ooooo 


Figure 3.5 

Dependency of the error of the calf angle on the separation between accelerometers 
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SOFTWARE 


4.1 INTRODUCTION 

The task achieved by the data acquisition program in C + + is to record synchronously 
the data from accelerometers and encoders, then to evaluate the angular positions, angular 
velocities and angular accelerations of the crank and pedals during exercise. A second program 
in MATLAB uses the data from the C + + program to determine the angular orientations, 
angular velocities, angular accelerations of the calf and thigh sections, and the linear 
accelerations of the knee joint and the hip joint. 

4.2 THE DATA ACQUISITION PROGRAM IN C+ + 

During the development of the software for data acquisition, five programs were 
created: SCANOPF.CPP, SCANOPF2.CPP, SCANOPB.CPP, SCANOPB2.CPP and 
SCANSTAR.CPP. Only SCANOPB.CPP determines angular positions, angular velocities and 
angular accelerations for the kinematics program in MATLAB. Each program served a purpose 
within the research. Figures 4.2— 4.6 show the flowcharts of the five programs. These programs 
differ from each other primarily in the way that they acquire and store data. The reason for 
creating the programs was to find a way to quickly digitize a large number of samples. Table 
4.1 shows the differences inside the data acquisition loops. The differences between 
SCANOPF.CPP, SCANOPB.CPP and SCANSTAR.CPP are in the ways that data is stored. The 
differences between SCANOPF.CPP and SCANOPF2.CPP, SCANOPB.CPP and 
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SCAN0PB2.CPP are that both SCANOPF2.CPP and SCANOPB2.CPP only acquire data from 
the A/D board but both SCANOPF.CPP and SCANOPB.CPP acquire data from the A/D board 
and the decoder board. The difference between the SCANSTAR.CPP and the other four is that 
the data buffer used in SCANSTAR.CPP is allocated by NI-DAQ but the buffers in the other 
four programs are allocated by DOS. The differences are explained in more detail later. 

4.2.1 Storage of Data 

The functions in the NI-DAQ library for scanning the channels are SCAN Op and 
SCAN_Start, and both can store data only into a buffer but not a file. The form of the data in 
the buffer is binary. This cannot be read directly from the screen, so it must be converted into 
ASCII text form (NI-DAQ Function Reference Manual). 

There are three possible methods to store data. The three methods are used by the five 
data acquisition programs, their differences are shown in Table 4.1 and their flowcharts appear 
in Figures 4.2-4. 6. The first method, used by SCANOPF.CPP and SCANOPF2.CPP, saves the 
acquired data from the buffer into a binary file immediately after each scan (i.e., save inside the 
data acquisition loop). Then the binary file is converted into a text file after the end of the 
entire data acquisition operation. The reason that the binary data in the buffer is not directly 
stored into a text file inside the data acquisition loop is to reduce the time inside the loop. The 
second method, used by SCANOPB.CPP and SCANOPB2.CPP, uses a buffer big enough for 
the entire data acquisition process, and then saves the data from the buffer into a text file. The 
difference between the first two methods is the size of the buffer used by the data acquisition 
operation. The first method requires a buffer only big enough for one scan but the second needs 
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Table 4.1 Differences inside data acquisition loop 


^^Differences 

Acquires 
data from 
the 

decoder 

board 

Buffer 


Files \ 

Size 

Allocated 

by 

A file to 
store data 

SCANOPF.CPP 

Yes 



Only big enough for one scan 

DOS 

Yes 

SCANOPF2.CPP 

No 

SCANOPB.CPP 

Yes 

Big enough to store the whole 

DOS 

No 

SCANOPB2.CPP 

No 

data of many scans 

SCANSTAR.CPP 

Yes 

Big enough to store the whole 
data of many scans 

NI-DAQ 

functions 

No 

1 


a buffer big enough for the entire data set. The third method is used by SCANSTAR.CPP and 
uses a buffer created by a NI-DAQ function, whose size can exceed 64k. 

The second method encounters the fact that Turbo C+ + (Version 3.0, Borland 
International Inc., Scotts Valley, CA) normally limits the size of all static data to 64K. The 
Huge Model sets aside that limit (Turbo C+ + User’s Guide). Yet the library functions offered 
by NI-DAQ are for the Large Model so they cannot work in the Huge Model. The first method 
avoids the buffersize limitation at the expense of reducing the rate of data acquisition, because 
it requires extra time in the data acquisition loop to store data into a file. The shortcoming of 
the third method is that the buffer is not available for synchronous data acquisition operations; 
thus, only SCAN_Start can use the third method. Therefore, only SCANOPB.CPP determines 
angular positions, angular velocities and angular accelerations, which are later used by the 
program for kinematics in MATLAB. 
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4.2.2 Synchronization with encoders 

Using the NI-DAQ function SCAN_Op, synchronous, multiple-channel scanned data 
acquisition operations are performed (NI-DAQ Function Reference Manual). Using SCANOp, 
synchronization is achieved by an internal on-board counter that produces trigger pulses. Thus 
A/D conversions do not begin until a trigger pulse is applied. This is controlled by two 
functions in NI-DAQ-DAQ_Config and CTR_Square. The former concerns the configuration 
information for subsequent data acquisition operations (NI-DAQ Function Reference Manual) 
and the latter programs a counter to generate a continuous square wave output of specified duty 
cycle and frequency (NI-DAQ Function Reference Manual). Therefore, the rate of trigger pulse 
of the counter defines the rate of operation of SCAN_Op. Each operation of SCAN Op is 
immediately followed by a reading of the decoder board. 
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Figure 4. 1 
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Figure 4.2 
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Figure 4.3 
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Figure 4.4 
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END ; 

FLOWCHART OF SACNSTAR.CPP 


Figure 4.5 
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4.2.3 Checking the D/A System 

From the algorithms shown in Figure 4.2 and Figure 4.4, the time difference between 
the starting of the data acquisition operation of the A/D board and the subsequent waiting of 
SCAN Op for the trigger pulse must be less than the period of the trigger pulse; otherwise, 
some trigger pulses may be missed. Whether trigger pulses are missed can be checked by 
inspecting a digitized signal whose original form was known already. If digitization does not 
affect the integrity of the signal, no trigger pulse is missed and the period is greater than the 
time difference. For example, given a sinusoidal signal at 5 Hz and the rate of the counter is 
50 Hz, the digitized sinusoidal signal should clearly reveal that ten data points are shown in a 
period if no trigger pulse is missed; otherwise, some trigger pulses are lost. Through this method 
the maximum rates of the entire data acquisition operation corresponding to different methods 
of storage were decided. The maximum rate of data acquisition in SCANOPF.CPP is 57Hz, 
the one of SCANOPB.CPP is 65Hz, and the one of SCANSTAR.CPP is above 1.2 kHz. 

4.2.4 Discussion 

The difference of the maximum rate of data acquisition between SCANSTAR.CPP and 
the other two is remarkable. It is necessary to analyze the essential difference between them to 
understand the cause of the big drop of the data acquisition rate. 

The main difference in the data acquisition algorithm of SCANSTAR.CPP and the other 
two is synchronization. SCANSATR.CPP does not read the decoder board. It only reads the DA 
board and was created as a research tool only. To synchronize the DA board and the decoder 



38 


board, SCANOPF.CPP and SCANOPB.CPP use a loop to include the data acquisition operation 
of the decoder board. Both the loop and the inclusion of the decoder board’s data acquisition 
operation reduce the data acquisition rate. If the main cause of the time delay was from the 
data acquisition operation of the decoder board, the big drop should be reduced by excluding the 
decoder board’s data acquisition operation from the data acquisition loop. If the exclusion of the 
decoder board’s data acquisition operation cannot reduce the big drop, the cause is in the call 
to SCANOp. 

After the decoder board’s data acquisition operation is excluded from the loops of 
SCANOPF.CPP and SCANOPB.CPP, the two programs become SCANOPF2.CPP and 
SCANOPB2.CPP, respectively. Figure 4.2 and Figure 4.4 show the flowcharts. The maximum 
rates of the two programs are 63Hz and 67Hz. This reveals that the main cause of the big drop 
of the data acquisition is due to the call to SCAN Op. 

A limiting factor in the system is that the most A/D conversions performed by any data 
acquisition function in NI-DAQ is 65535. This limit prevents long data acquisition time at high 
rate. The programs that save data to a file after every scan can acquire data as long as there is 
space in the hard disk, but the rate may be too slow for some applications. 

4.3 THE MATLAB PROGRAM FOR KINEMATICS 

The program, KINET.M, uses the parameters 0,,0 2 ,a,,o! 2 ,c<Ji,a) 2 and accelerometer values 
from an external file made by SCANOPB.CPP to determine the angular positions, angular 
velocities, angular accelerations and linear accelerations of the calf and thigh sections during 
pedaling. Its flowchart is shown in Figure 4.7. Section 3.3, Kinematic analysis, has described 
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the theoretical background of the program. 


(|tart) 

Load 0,w, a of crank and pedals, and the readings of accelerometers 
from an external file produced by a data acquisition program in C++ 


Determine 0 a 2 using the 0, w and a 

i 

Determine the co and a using the readings of the accelerometers on calf 

3 3 


Determine 3 339 


Determine 3 g 


Determine 6^ using °a 2 and a 2 g 


Determine usin 9 ®03 and a j 


Repeat the preceding procedure to determine oo 4 , a 4 , and 0 M 


END 


FLOWCHART OF KINET.M 


Figure 4.6 The flowchart of KINET.M 


4.4 SIMULATION 

Using as input 0, ,^ 2 ,a 1 ,a 2 ,oo 1 ,oo 2 and the accelerometers’ measurements, KINET.M 
determines <23,0*4,003,004 and the rotation angles between the moving frames and the inertial 
frame. A simulation was done to make sure that KINET.M works appropriately. 

Many sets of test values were used as inputs to a test program, TEST.M. Every set of 
inputs contained 0 U d 2 , a t , a 2 , co,, co 2 , 0 o 3 , 0^, a 3 , o* 4 , oo 3 and oo 4 . TEST.M produced consistent 
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accelerometer values according to the given inputs. These values were then used as inputs to 
KINET.M to determine 0 O3 , 0^, a 3 , a 4 , w 3 and oj 4 . Using the same 0,, 0 2 , a„ a 2 , w, and o> 2 and 
the accelerometer values from TEST.M, KINET.M produced identical 0 O3 , 004, o; 3 , a 4 , u 3 and 
w 4 to the inputs to TEST.M. 

In TEST.M °a 2 was first determined using Equation 3.1, which was then used to 
determine 3 a 2 using Equation 4.1. 

fcos(0 O 3) -sinCOos)] -1 V (4 ^ 

\ sin(0 O3 )' cos(0 O3 )J \ y 

The following two equations are used to determine 3 a 3 , 


3 «Zx + 3, *32y a 3 + 3r 32x W 3 2 


3 «3, = 3fl 2, " 3r 32x a 3 - 3r 32y W 3 2 (4 ’ 3) 


where 3 r 32 stands for the position vector of the foot joint w.r.t Frame 3. 

Finally the theoretical values of the accelerometers attached to the calf were determined 


as follows. 



(4.4) 



(4.5) 


where 3 r cl and 3 r c2 represent the position vectors of accelerometers attached to the calf with 
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respect to Frame 3. 

Since gravity exerts a -g acceleration on accelerometers, the theoretical accelerometer 
values were converted to their respective measurements as follows (Equations 4.6 and 4.7). 




gsin(% 3 ) 

gcos(Q Q3 ) 


(4.6) 




gsw(0 O3 ) 

SC0s(0 O3 ) 


(4.7) 


For the thigh section, °a 3 was first obtained using Equation 4.8. 


K [ c °s(e 03 ) -sm(0 O3 )l ^ (4 g) 

°a 3 J [sin(0 O3 ) cos(0 O3 ) J [ \ y 
This procedure was repeated for the calf section to produce the measurements of the 
accelerometers attached to the thigh section, a,, and a^. 

The inputs to TEST.M, the consistent accelerometer values produced by TEST.M. and 
the output of KINET.M are shown in Appendix A.1-A.3, respectively. The first six rows in 
Appendix A.1-A.3 are all the same and represent the designed 0,, <*>,, a,, 0 2 , a> 2 , and os, 
respectively. The last four rows in Appendix A.l are the designed 0 O3 , co 3 , a 3 , w 4 and a 4 . 

Accelerometer values produced by TEST.M are the last eight rows representing a clx , a clv , a c3x , 
a c3y , a, lx , a tl y, a^ and a^y, respectively. The results of simulation shown in Appendix A. 3 
reveal that the kinematics program works perfectly. 



CONCLUSIONS 


A two dimensional analytical model to determine the kinematic parameters of a subject 
bicycling was developed based on a previously specified measurement system. It was 
implemented in a MATLAB program, and then proved to work by simulation. This system and 
method may be used to determine the angular velocity and angular acceleration, linear 
acceleration, position and orientation of a moving object, knowing the acceleration and position 
of a known point in the object with respect to the base frame. The precision error of 
accelerometer values results in an error of the rotation angle between the inertial frame and the 
moving frame attached to the lower limbs. A theoretical analysis of the resultant error offers an 
estimate and methods to reduce the error. Reduction of the error can be achieved by 
increasing the distance between the accelerometers. The results of the numerical error 
calculations based on the simplified five-bar linkage model are shown in Figures 3.4 and 3.5. 
These figures show that the errors did decrease with larger distance between the accelerometers. 
The maximum errors corresponding to 5% accelerometer precision error were determined to be 
about 6 and 8.25 degrees for different separations between accelerometers; and the maximum 
errors corresponding to 3% accelerometer precision error were determined to be about 3.6 and 
5 degrees. 

The A/D board and the decoder board have been proven to work well. The 
accelerometers are not yet fully usable, and that is why an experimental verification is outside 
the objective of this thesis. The encoders have been tested and work well. The five D/A 
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programs have also been proven to work well. The fastest rates of data acquisition of the 
programs SCANOPF.CPP, SCANOPF2.CPP, SCANOPB.CPP, SCANOPB2.CPP and 
SCANSTAR.CPP are 57 Hz, 63 Hz, 65 Hz, 67 Hz and above 1.2 kHz, respectively. A 
complete experimental verification of the proposed system cannot be realized until the 
hardware is ready. 

To quantify bone loading during exercise, forces exerted by a particular muscle, or 
muscle group, and forces and torques at joints need to be determined. The measurement system 
can be extended to include load cells to measure contact forces, and EMG probes to provide 
qualitative information about muscle action. Using these additional sensors, and optimization 
techniques, it is possible to determine musculoskeletal loading. 
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RECOMMENDATIONS 


Because of the need to synchronize with the decoder board’s data acquisition operation, 
a large drop of the rate of data acquisition occurs. The drop should be greatly reduced if the 
A/D board can directly acquire the encoders’ data. Usually, encoder data includes two chains 
of TTL pulses which are decoded by the decoder board. Therefore, a program needs to be 
developed to interpret encoder data that would bq digitized using the A/D board (Philippe 
Coiffet, and Michel Chirouze, 1982). 

The model for determining kinematic parameters in this research only works for two 
dimensions. This model should be generalized to the three dimensional case. The generalization 
is left for future work. 
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APPENDIX A. 1 INPUTS TO TEST.M 



input = 


Columns 1 through 7 


10 .0000 

10.0000 

10.0000 

10 . 0000 

10 . 0000 

10 . 0000 

30.0000 

2 .0000 

2 .0000 

2.0000 

2 .0000 

2.0000 

2.0000 

2 .0000 

o 

0 

0 

0 

0 

0 

0 

30 . 0000 

30 .0000 

30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

3 .0000 

2.0000 

1.0000 

3.0000 

2 . 0000 

1.0000 

3 . 0000 

1 .5000 

0 

0 

1.5000 

0 

0 

1.5000 

4 .0000 

1.0000 

2.0000 

4 . 0000 

1 . 0000 

2 .0000 

4 . 0000 

2 .5000 

3.0000 

7.7330 

2.5000 

3.0000 

7.7330 

2.5000 

45 .0000 

50.0000 

55.0000 

60.0000 

65 . 0000 

70 . 0000 

75.0000 

2.0000 

3.0000 

4 . 0000 

2.0000 

3 . 0000 

4.0000 

2 .0000 

1.0000 

7 .0000 

9.0000 

1.0000 

7.0000 

9.0000 

1 .0000 

135.0000 

140.0000 

145.0000 

150.0000 

155.0000 

160.0000 

165.0000 

Columns 8 

through 14 





30.0000 

30.0000 

30 . 0000 

30.0000 

30.0000 

30.0000 

30.0000 

2.0000 

2.0000 

2.0000 

2.0000 

2 . 0000 

2.0000 

2 . 0000 

0 

0 

0 

0 

0 

0 

0 

30 .0000 

30 . 0000 

30.0000 

30.0000 

30 . 0000 

30 . 0000 

30.0000 

2 . 0000 

1 . 0000 

3 . 0000 

2 .0000 

1.0000 

3.0000 

2.0000 

0 

0 

1.5000 

0 

0 

1.5000 

0 

1 . 0000 

2.0000 

4.0000 

1.0000 

2 .0000 

4.0000 

1.0000 

3 .0000 

7.7330 

2.5000 

3.0000 

7.7330 

2.5000 

3 .0000 

80.0000 

85.0000 

90.0000 

95.0000 

100 . 0000 

105.0000 

110 . 0000 

3 .0000 

4.0000 

2.0000 

3.0000 

4.0000 

2.0000 

3.0000 

7.0000 

9.0000 

1.0000 

7 . 0000 ' 

9.0000 

1 .0000 

7.0000 

170.0000 

175.0000 

180 . 0000 

185.0000 

190.0000 

195.0000 

200 .0000 

Columns 15 through 

19 





30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

• 


2 . 0000 

2.0000 

2.0000 

2.0000 

2.0000 



0 

0 

0 

0 

0 



30.0000 

30 . 0000 

30.0000 

30.0000 

30.0000 



1 .0000 

3.0000 

2.0000 

1.0000 

1.0000 



0 

1.5000 

0 

0 

0 



2.0000 

4.0000 

1.0000 

2.0000 

2 .0000 



7.7330 

2 .5000 

3 . 0000 

7.7330 

7 . 7330 



115 .0000 

120.0000 

125.0000 

130.0000 

135 . 0000 



4.0000 

2.0000 

3.0000 

4.0000 

4.0000 



9.0000 

1.0000 

7.0000 

9.0000 

9.0000 



205.0000 

210.0000 

215 . 0000 

220.0000 

225.0000 





APPENDIX A. 2 ACCELEROMETER VALUES PRODUCED BY TEST.M 



accelerometer 


1.0e+003 * 


Columns 1 through 7 


0.0002 

0.0002 

0.0002 

0.0002 

0.0002 

0.0002 

0.0005 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0 

0 

0 

0 

0 

0 

0 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0030 

0.0020 

0.0010 

0.0030 

0.0020 

0.0010 

0.0030 

0.0015 

0 

0 

0.0015 

0 

0 

0.0015 

-1.6468 

-1.1093 

-1.1135 

-1.7509 

-1.2010 

-1.1891 

-1.8695 

-0.4026 

-0.4033 

-0.1762 

-0.0784 

-0 . 1031 

0.1056 

0.1191 

-1.3268 

-1.0893 

-1.0335 

-1.4309 

-1.1810 

-1.1091 

-1.5495 

-0.4526 

-0.4633 

-0.3308 

-0.1284 

-0.1631 

-0.0491 

0.0691 

-0.5726 

-0.6133 

-0.5538 

-0.2484 

-0.3131 

-0.2721 

-0.0509 

1 . 8 653 

1.3893 

1.5422 

1.9734 

1.4810 

1.6178 

2.0920 

-0.4926 

-0.4333 

-0.2338 

-0.1684 

-0.1331 

0.0479 

0.0291 

1.8493 

1.2493 

1.3622 

1.9534 

1.3410 

1.4378 

2.0720 

Columns 8 

through 14 






0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0.0020 

0 

0 

0 

0 

0 

0 

0 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0 .0005 

0.0005 

0.0020 

0.0010 

0.0030 

0.0020 

0 .0010 

0.0030 

0.0020 

0 

0 

0.0015 

0 

0 

0.0015 

0 

-1.2956 

-1.2603 

-1.8310 

-1.24 >0 

-1.1977 

-1.6996 

-1.1145 

0. 1142 

0.3383 

0.4832 

0.4450 

0.6405 

0.8249 

0.7519 

-1.2756 

-1.1803 

-1.5110 

-1.2270 

-1.1177 

-1.3796 

-1.0945 

0.0542 

0.1837 

0.4332 

0.3850 

0.4858 

0.7749 

0 . 6919 

-0.0958 

-0.0394 

0.3132 

0.2350 

0.2628 

0.6549 

0.5419 

1.5756 

1.6889 

2.0535 

1.5270 

1 . 6264 

1 9 2 2 1 

1.3945 

0.0842 

0.2806 

0.3932 

0.4150 

0.5828 

0.7349 

0 . 7219 

1 . 4356 

1.5089 

2.0335 

1.3870 

1.4464 

1.9021 

1.2545 

Columns 15 

through 19 






0.0005 

0.0005 

0.0005 

0.0005 

0 . 0005 



0.0020 

0.0020 

0.0020 

0.0020 

0.0020 



0 

0 

0 

0 

0 



0.0005 

0.0005 

0.0005 

0.0005 

0 .0005 



0.0010 

0.0030 

0.0020 

0.0010 

0.0010 



0 

0.0015 

0 

0 

0 



-1.0592 

-1.4842 

-0.9070 

-0.8539 

-0.7730 



0.9162 

1.1210 

1.0141 

1.1467 

1.2106 



-0.9792 

-1.1642 

-0.8870 

-0.7739 

-0.6930 



0.7615 

1.0710 

0 . 9541 

0.9920 

1.0559 



0.5385 

0.9510 

0.8041 

0.7690 

0.8329 



1.4878 

1.7067 

1.1870 

1.2826 

1.2017 



0.8585 

1.0310 

0.9841 

1.0890 

1.1529 



1.3078 

1.6867 

1.0470 

1.1026 

1.0217 





APPENDIX A. 3 THE RESULTS OF SIMULATION 



simulation = 

Columns 1 

10.0000 

2.0000 

0 

30.0000 

3.0000 

1.5000 
4 .0000 

2.5000 

45.0000 

2.0000 

1.0000 

135 .0000 

Columns 8 

30.0000 

2.0000 

0 

30.0000 

2.0000 

0 

1.0000 

3.0000 

80.0000 

3.0000 

7.0000 

170.0000 

Columns 15 

30.0000 

2.0000 
0 

30.0000 

1.0000 

0 

2.0000 

7.7330 

115.0000 

4.0000 

5 .0000 

205 .0000 


through 7 


10.0000 

10.0000 

10.0000 

10.0000 

10.0000 

30.0000 

2.0000 

2.0000 

2.0000 

2.0000 

2.0000 

2.0000 

0 

0 

0 

0 

0 

0 

30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

2.0000 

1 . 0000 

3.0000 

2.0000 

1.0000 

3.0000 

0 

0 

1.5000 

0 

0 

1.5000 

1.0000 

2.0000 

4.0000 

1.0000 

2.0000 

4.0000 

3.0000 

7.7330 

2.5000 

3.0000 

7.7330 

2.5000 

50.0000 

55.0000 

60.0000 

65.0000 

70.0000 

75.0000 

3 . 0000 

4.0000 

2.0000 

3.0000 

4.0000 

2.0000 

7 . 0000 

9.0000 

1.0000 

7 .0000 

9.0000 

1.0000 

140.0000 

145.0000 

150.0000 

155.0000 

160.0000 

165.0000 

through ! 

14 





30 . 0000 

30.0000 

30.0000 

30 . 0000 

30.0000 

30.0000 

2 . 0000 

2.0000 

2.0000 

2.0000 

2.0000 

2.0000 

0 

0 

0 

0 

0 

0 

30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

30.0000 

1.0000 

3.0000 

2.0000 

1.0000 

3.0000 

2.0000 

0 

1.5000 

0 

0 

1.5000 

0 

2.0000 

4.0000 

1.0000 

2 . 0000 

4.0000 

1.0000 

7 . 7330 

2.5000 

3.0000 

7.7330 

2.5000 

3.0000 

85.0000 

90.0000 

95.0000 

100.0000 

105.0000 

110.0000 

4.0000 

2.0000 

3.0000 

4.0000 

2.0000 

3.0000 

9.0000 

1.0000 

7 . 0000 * 

~ 9.0000 

1.0000 

f 7.0000 

175.0000 

180.0000 

185.0000 

190.0000 

195.0000 

200.0000 

through 

19 





30.0000 

30.0000 

30.0000 

30.0000 

* 


2.0000 

2.0000 

2.0000 

2.0000 



0 

0 

0 

0 



30 . 0000 

30 . 0000 

30.0000 

30.0000 



3.0000 

2 . 0000 

1.0000 

1.0000 



1.5000 

0 

0 

0 



4.0000 

1 . 0000 

2.0000 

2 .0000 



2 .5000 

3 . 0000 

7.7330 

7.7330 



120.0000 

125.0000 

130.0000 

135 .0000 



2.0000 

3.0000 

4.0000 

4 .0000 



1.0000 

7.0000 

9.0000 

9.0000 



210.0000 

215.0000 

220.0000 

225.0000 





APPENDIX B.l KINET.M 



^ ************** ****** ******** *********************************** 
% set the parameters in file and read from it to determine variable 

or *************************************************************** 


function x=kinet(kar,b); 

% b is a scale factor 
% define some constant variable 

11= 15*b ; % crank length 

12= 10*b ; % the length between foot joint and 

% pedal 


13= 40*b ; % calf 

14= 40*b ; % thigh 


% % % % % % % % % % for left pedal 

rcl=[ -10,-5]’*b ; % in frame 3 

rc2=[ -30,-5]’*b ; % in frame 3 

rtl = [ -10,-5]’*b ; % in frame 4 

rt2 = [ -30,-5]’ *b ; % in frame 4 

r32=[ -13 0]’; 
r43=[ -14 0]’; 


% for right pedal 

% rcl=[ -2,-2]’*b ; % in frame 3 

% rc2=[ -4,-2] ’*b ; % in frame 3 

% rtl=[ -2,-2]’*b ; % in frame 4 

% rt2 = [ -4,-2] ’*b ; % in frame 4 

% r32=[ -13 0]’ ; 

% r43=[ -14 0]’ ; 


g=980*b; 



% get paremetesr from data acquisition program 


[m,n] = size(kar); 

crank(l,:)=kar(l,:); 
w_c(l,:)=kar(2,:); 
arfa_c( 1 , :) =kar(3 , : ) ; 

%%%%%%%%%%%%%%%%%% loop for kinetmatics %%%%%%%%%%%%% 


for j = 1:2 % 2 for right pedal, position vectors need to be redefined if not the same as 

% the left 

* ^ ► 

if j = = l, %%%%%%%%%%for left pedal 

pedal(l,:)=kar(4,:); 

w_p(l,:)=kar(5,:); 

arfa_p(l,:)=kar(6,:); 

acl(l:2,:)=kar(7:8,:); 

atl(l:2,:)=kar(ll:12,:); 

ac2(l:2,:)=kar(9: 10,:); 

at2(l:2,:)=kar(13: 14,:); 


else %%%%%%% %%%%for right pedal 

pedal(l,:)=kar(15,:); 
for i = l:n 

crank( 1 , i) = crank( 1 , i) + p i ; 
end 

w_p(l,:)=kar(16,:); 
arfa_p(l,:)=kar(17,:); 
acl(l:2,:)=kar(18: 19,:); 
atl(l:2,:) = kar(22:23,:); 
ac2(l:2,:)=kar(20:21,:); 
at2(l:2,:)=kar(24:25,:); 

end 

rcl2=rc2-rcl; 

rtl2=rt2-rtl; 



for i = l:n 


acl2=ac2(:,i)-acl(:,i); 

atl2=at2(:,i)-atl(:,i); 


Sl=sin(crank(l,i)); 

Cl=cos(crank(l,i)); 


% determine a02 in frame 0 

S 12 = sin(crank( 1 , i) 4- pedal( 1 , i)) ; 

C 1 2 = cos(crank( 1 , i) + pedal ( 1 , i)) ; 
t3=[-ll*Sl -12*S12; 11*C1 12*C12]; 

t4=[ arfa_c(l,i) (arfa_c(l,i)+arfa_p(l,i))]’; 

t5=[ 11*C1 12*02; 

11*S1 12*S12]; 

t6=[ (w_c(l,i)"2) (w_c(l,i)+w_p(l,i)) A 2]’; 

a02(l:2,i)=t3*t4-t5*t6; % in frame 0 


% determine a03 in frame3 

tl =inv(arr(rcl2)) :i 'acl2; 
w_ca( 1 , i) = sqrt(t 1 (2,1)); 
arfa_ca(l,i)=tl(l,l); 
t2=arr(rcl)*tl; 

a03_3(l:2,i)=-t2+acl(:,i); % in frame 3 


% determine rotation angle between frame 0 and frame 3 



a02 3=a03 3(:,i)+arr(r32)*tl; % in frame 3 and without compensation of 

% gravity 

% a02 3 should be a02_3+com(g,theat03) 


t8(i)=(a02(l,i)*a02_3(U)+a02(2,i)^a02_3(2,l)-a02_3(2,l)*g)/((a02_3(l,l) A 2) + (a02_3(2,l) 

" 2 »; 

t9=(a02(2,i)*a02_3(l, l)-a02(l ,i)*a02_3(2, l)-a02_3(l , l)*g)/((a02_3(l , 1)^2) +(a02_3(2, 1) A 2)) 


theta03( 1 , i) =acos(t8(i)) ; 


a03 3(1:2, i)=a03_3(l:2,i)+commm(g,theta03(:,i)); % COMPENSATION WITH 

% GRAVITY; 


a03(l:2,i)=(rot2d(theta03(:,i)))*a03_3(:,i); % in frame 0; 

%%%%%%% % determine a04 in frame4 %%%%%%%%%%%%%%%%%%%% 


tl = inv(arr(rtl2))*atl2; 
w_t(l,i)=sqrt(t 1(2,1)); 
arfa_t(l,i)=tl(l,l); 
t2=arr(rtl)*tl; 

a04_4(:,i)=-t2+atl(:,i); % in frame 4 

% determine rotation angle between frame 0 and frame 4 
a 03_4 = a04_4 (:,i)+arr(r43)*tl; % in frame 3 and without compensation of 

% gravity 

% it should be a03_4+com(g,theta04) 


theta04(l,i) =acos((a03(l ,i)*a03_4(l , 1) +a03(2,i)*a03_4(2, l)-a03_4(2, l)*g)/((a03_4(l , 1)"2) 
+(a03_4(2, 1)^2))); 


t9 = (a03(2,i)*a03_4(l,l)-a03(l,i)*a03_4(2,l)-a03_4(l,l)*g)/((a03_4(l,l)^2) + (a03_4(2,l)^2)) 



if t9 < =0 


theta04(l,i)=pi+asin(-t9); 

end 


a04_4(:,i)=a04_4(:,i)+commm(g,theta04(:,i)); 
a04(: ,i) =(rot2d(theta04(: ,i)))*a04_4(: ,0; 
end % end of operation loop 


% COMPENSATION WITH 

% GRAVITY 


if j = = l 

a03_l(:,:)=a03(:,:); 
a04_l(:,:)=a04(:,:); 
theta03_l(l , :) = theta03( 1 , :) : * 1 80/pi; 
theta04_I( 1 ,:)= theta04( 1 ,:) * 1 80/pi; 
arfa_cal ( 1 , : ) = ar fa_ca( 1 , : ) ; 
w_cal( 1 , :) = w_ca( I , : ) ; 
arfatl ( 1 , : ) = ar fa_t( 1 , : ) ; 
w_tl(l,:)=w_t(l,:); 


else 

a03_r(:,:)=a03(:,:) 

a04_r(:,:)=a04(:,:) 

theta03_r ( l,:)=theta03(l,:)/pi 5|t 180; 

theta04_r ( 1 , : ) = theta04 (l,:)/pi*180; 

arfa_car( 1 , :) = arfa_ca( 1 , : ) ; 

w_car(l,:)=w_ca(l,:); 

arfatr ( 1 , : ) = ar fa_t( 1 , : ) ; 

w_tr(l,:)=w_t(l,:); 

end % end of if 


end % end for j 



x(l,:)=w_cal(l,:); 
x(2 , : ) =arfa_cal( 1 , : ) ; 
x(3 , : ) =■ theta03 _1 ( 1 , : ) 
x(4,:)=w_tl(l,:); 
x(5,:)=arfa_tl(l,:); 
x(6,:)=theta04_l(l,:) 
x(7:8,:)=a03_l(:,:); 
x(9:10,:)=a04_l(:,:); 

x(ll,:)=w_car(l,:); 

x(12,:)=arfa_car(l,:) 

x( 1 3 , : ) = theta03_r ( 1 , 

x(14,:)=w_tr(l,:); 

x(15,:)=arfa_tr(l,:); 

x(16,:)=theta04_r(l, 

x(17:18,:)=a03_r(:,: 

x(19:20,:)=a04_r(:,: 


* 180/pi; 
80/pi; 

:’)* 180/pi; 
:)* 1 80/pi; 



APPENDIX B.2 TEST.M 



function l=test(kar,b); 


% define some constant variable 

11= 15*b ; % crank length 

12= 10*b ; % the length between foot joint and 

% pedal 

13= 40*b ; % calf 

14= 40*b ; % thigh 

%%%%%%%%%% for left pedal 

rcl=[ -10,-5]’*b ; % in frame 3 

rc2 = [ -30,-5]’ *b ; % in frame 3 

rtl=[ -10,-5]’*b ; % in frame 4 

rt2=[ -30,-5] ’*b ; % in frame 4 

r32=[ -13 0]’; 
r43 = [ -14 0]’; 

g=980*b; 

% get paremetesr from data acquisition program 


[m,n]=size(kar); 

crank(l,:)=kar(l,:)*pi/180; 

w_c(l,:)=kar(2,:); 

arfa_c( 1 , :) = kar(3 , :) ; 

pedal ( 1 , : ) = kar (4,:)*pi/180; 

w_p(l,:)=kar(5,:); 

arfa_p(l,:)=kar(6,:); 

w_ca(l,:)=kar(7,:); 

theta03(l,:)=kar(9,:)*pi/180; 

ar fa_ca( 1 , : ) = kar (8 , : ) ; 

w_t(l,:)=kar(10,:); 

arfa_t(l,:)=kar(l 1,:); 

theta04( 1 , : ) = kar ( 1 2 , : ) *p i/ 1 80 ; 


for i = l:n 



% determined aOl in frame 0 
SI =sin(crank(l,i)); 

Cl =cos(crank(l,i)); 


% determine a02 in frame 0 

S 12 =sin(crank( 1 , i) +pedal( 1 , i)); 

C 1 2 = cos(crank( 1 , i) +pedal ( 1 , i)) ; 
t3 = [ -11*S1 -12*S12; 11*0 12*C12]; 

t4 = [ arfa_c(l,i) (arfa_c(l,i)+arfa_p(l,i))]’; 

t5 = [ 11*0 12*C12; 

11*S1 12*S12]; 

t6 = [ (w_c(l,i) A 2) (w_c(l,i)+w_p(l,i)) A 2]’; 

a02(l:2,i)=t3*t4-t5*t6; % in frame 0 

k = [arfa_ca( 1 , i ) w_ca( 1 , i) A 2 ] ’ ; 

a02_3(l :2, i) = inv(rot2d(theta03( 1 , i)))*a02(: , i); 

a03_3( 1 :2, i) =a02_3( 1 :2, i)-arr(r32)*k ; 

ac 1 ( : , i) = a03_3 ( : , i) + arr (rc 1 ) *k-co m m m(g , theta03 ( : , i)) ; 

ac2(:,i)=a03_3(:,i)+arr(rc2)*k-commm(g,theta03(:,i)); 

a03(: , i) = rot2d(theta03( 1 , i)) *a03_3 ( : , i) ; 

k = [arfa_t( 1 , i) w_t(l,i) A 2]’; 

a03_4 ( 1 : 2 , i) = i nv (rot2 d (the ta04 ( 1 , i))) *a03 ( : , i) ; 

a04_4( 1 : 2 , i) = a03_4( 1:2, i)-arr(r43) *k; 

at 1 (: , i) = a04_4( : , i) + arr(rt 1 ) 5k k-commm(g , theta04(: , i)) ; 

at2(:,i)=a04_4(:,i)+arr(rt2) 5tc k-commm(g,theta04(:,i)); 

end 


l(l,:)=crank(l,:); 
I(2,:)=w_c(l,:); 
l(3,:)=arfa_c(l,:); 
l(4,:)=pedal(l,:); 
1(5,:)= w_p(l,:); 
l(6,:)=arfa_p(l,:); 
1(7:8,:)= ac 1(1:2,:); 
1(9:10, :)=ac2(l:2,:); 
1(11:12, :)=atl(l:2,:); 
1(13: 14,:)= at2 (1:2,:); 



